

## load packages and data
pacman::p_load(estimatr, tidyverse, texreg, dotwhisker, data.table)

load("Data/Replication/Miller/Miller_Replication_Ready.rda")




## prepare depvar, indepvar combinations
depvars <- c("infant", "under5", "life", "literacy_original", "schooling2", "lit_equ", "school_equ")

rxdepvars <- c("rx_infant", "rx_under5", "rx_life", "rx_literacy_original", "rx_schooling2", "rx_lit_equ", "rx_school_equ")


indepvars <- c("ea" , "dem",
               "RegType_lied_EA", "RegType_lied_DEM",
               "RegType_RoW_EA", "RegType_RoW_DEM",
               "Politytype_Anocracy", "Politytype_Democracy",
               "status_fh_PF", "status_fh_F",
               "HTW_RegType_MP_Autocracy", "HTW_RegType_Democracy",
               "RegType_magaloni_EA", "RegType_magaloni_DEM",
               "AnckarRegtype_MP_Autocracy", "AnckarRegtype_Democracy")


controls <- c("aid", "loggdp",  "gdp_grow", "gini2",  "resourcesdep_hm", 
              "communist", "urban_cow", "ELF", "lnpop", "violence_domestic")


## run all models
a <- seq(1,15 , by = 2)
models <- list()


for (i in 1:length(depvars)) {
  
  outcome <- depvars[i]
  rx <- rxdepvars[i]
  
  for (j in 1:8) {
    print(paste0("Running Model " , outcome,"_model_",j, ", you bastard!"))
    
    indepvar <- indepvars[c(a[j],a[j]+1)]
    
    formula <- as.formula(paste(outcome, "~",  paste(indepvar, collapse = "+"), "+", rx, "+", paste(controls, collapse = "+")))
    
    
    models[[paste0(outcome,"_model_", j)]] <- lm_robust(formula, data = df, fixed_effects = year, clusters = ccode, se_type = "stata")
    
  }
}


submodels <- list()

for (i in 1:length(models)){
  submodels[[paste0("sub_", models[i])]] <- tidy(models[[i]]) %>%
    filter(term %in% c("ea", "RegType_lied_EA", "RegType_RoW_EA", "Politytype_Anocracy", "status_fh_PF", "AnckarRegtype_MP_Autocracy", "HTW_RegType_MP_Autocracy", "RegType_magaloni_EA")) %>% 
    mutate(model = names(models)[i])
}

plotmodels <- submodels %>% 
  rbindlist %>% 
  as.data.frame
rm(submodels)


plotmodels$term2 <- plotmodels$term
plotmodels$model2 <- plotmodels$model
plotmodels$term <- plotmodels$outcome
plotmodels$model <- plotmodels$term2

unique(plotmodels$term)
plotmodels <- plotmodels%>%mutate(model2 = case_when(
  model == "ea" ~ "EA Original",
  model == "RegType_lied_EA" ~ "LIED",
  model == "RegType_RoW_EA" ~ "RoW",
  model == "Politytype_Anocracy" ~ "Polity",
  model == "status_fh_PF" ~ "FH",
  model == "AnckarRegtype_MP_Autocracy" ~ "CPR",
  model == "HTW_RegType_MP_Autocracy" ~ "ARD",
  model == "RegType_magaloni_EA" ~ "AoW"))
plotmodels$model <- plotmodels$model2


plotmodels <- plotmodels%>%mutate(term2 = case_when(
  term == "under5" ~ "Under-5 Mortality",
  term == "infant" ~ "Infant Mortality",
  term == "life" ~ "Life Expectancy",
  term == "literacy_original" ~ "Literacy",
  term == "schooling2" ~ "Schooling",
  term == "lit_equ" ~ "Literacy Equality",
  term == "school_equ" ~ "Schooling Equality"))

plotmodels$term <- plotmodels$term2



## New plot

plotmodels<- plotmodels %>% mutate ( conf.low_95 = estimate - 1.96 *  std.error, conf.high_95 = estimate + 1.96 * std.error,
                                     conf.low_90 = estimate - 1.64 *  std.error, conf.high_90 = estimate + 1.64 * std.error) %>% as_tibble()






### add regime history
indepvars <- c("ea_avg" , "dem_avg",
               "RegType_lied_EA_hist", "RegType_lied_DEM_hist",
               "RegType_RoW_EA_hist", "RegType_RoW_DEM_hist",
               "Politytype_Anocracy_hist", "Politytype_Democracy_hist",
               "status_fh_PF_hist", "status_fh_F_hist",
               "HTW_RegType_MP_Autocracy_hist", "HTW_RegType_Democracy_hist",
               "RegType_magaloni_EA_hist", "RegType_magaloni_DEM_hist",
               "AnckarRegtype_MP_Autocracy_hist", "AnckarRegtype_Democracy_hist"
)


depvars <- c("infant", "under5", "life", "literacy_original", "schooling2", "lit_equ", "school_equ")

rxdepvars <- c("rx_infant", "rx_under5", "rx_life", "rx_literacy_original", "rx_schooling2", "rx_lit_equ", "rx_school_equ")




controls <- c("aid", "loggdp",  "gdp_grow", "gini2",  "resourcesdep_hm", 
              "communist", "urban_cow", "ELF", "lnpop", "violence_domestic")



a <- seq(1,15 , by = 2)
models <- list()


for (i in 1:length(depvars)) {
  
  outcome <- depvars[i]
  rx <- rxdepvars[i]
  
  for (j in 1:8) {
    print(paste0("Running Model " , outcome,"_model_",j, ", you bastard!"))
    
    indepvar <- indepvars[c(a[j],a[j]+1)]
    
    formula <- as.formula(paste(outcome, "~",  paste(indepvar, collapse = "+"), "+", rx, "+", paste(controls, collapse = "+")))
    
    
    models[[paste0(outcome,"_model_", j)]] <- lm_robust(formula, data = df, fixed_effects = year, clusters = ccode, se_type = "stata")
    
  }
}

submodels <- list()

for (i in 1:length(models)){
  submodels[[paste0("sub_", models[i])]] <- tidy(models[[i]]) %>%
    filter(term %in% c("ea_avg", "RegType_lied_EA_hist", "RegType_RoW_EA_hist", "Politytype_Anocracy_hist", "status_fh_PF_hist", "AnckarRegtype_MP_Autocracy_hist", "HTW_RegType_MP_Autocracy_hist", "RegType_magaloni_EA_hist")) %>% 
    mutate(model = names(models)[i])
}

plotmodels2 <- submodels %>% 
  rbindlist %>% 
  as.data.frame
rm(submodels)


plotmodels2$term2 <- plotmodels2$term
plotmodels2$model2 <- plotmodels2$model
plotmodels2$term <- plotmodels2$outcome
plotmodels2$model <- plotmodels2$term2

unique(plotmodels2$term)
plotmodels2 <- plotmodels2%>%mutate(model2 = case_when(
  model == "ea_avg" ~ "EA Original",
  model == "RegType_lied_EA_hist" ~ "LIED",
  model == "RegType_RoW_EA_hist" ~ "RoW",
  model == "Politytype_Anocracy_hist" ~ "Polity",
  model == "status_fh_PF_hist" ~ "FH",
  model == "AnckarRegtype_MP_Autocracy_hist" ~ "CPR",
  model == "HTW_RegType_MP_Autocracy_hist" ~ "ARD",
  model == "RegType_magaloni_EA_hist" ~ "AoW"
))
plotmodels2$model <- plotmodels2$model2

plotmodels2 <- plotmodels2%>%mutate(term2 = case_when(
  term == "under5" ~ "Under-5 Mortality",
  term == "infant" ~ "Infant Mortality",
  term == "life" ~ "Life Expectancy",
  term == "literacy_original" ~ "Literacy",
  term == "schooling2" ~ "Schooling",
  term == "lit_equ" ~ "Literacy Equality",
  term == "school_equ" ~ "Schooling Equality"))

plotmodels2$term <- plotmodels2$term2

plotmodels2<- plotmodels2 %>% mutate ( conf.low_95 = estimate - 1.96 *  std.error, conf.high_95 = estimate + 1.96 * std.error,
                                     conf.low_90 = estimate - 1.64 *  std.error, conf.high_90 = estimate + 1.64 * std.error) %>% as_tibble()




plotmodels2$bla <- "Coef"
plotmodels2$type<- "Regime History"


plotmodels$bla <- "Coef"
plotmodels$type<- "Regime Dummy"



names(plotmodels)
names(plotmodels)

plotframe <- rbind(plotmodels, plotmodels2)




plotframe$term <- fct_rev(factor(plotframe$term, levels=unique(plotframe$term)))

plotframe$model <- fct_rev(factor(plotframe$model, levels=unique(plotframe$model)))

unique(plotframe$term)



### Make Figure K1
plotframe %>% 
  filter(term %in% c("Infant Mortality", "Life Expectancy", "Schooling", "Schooling Equality"))%>%
  ggplot(aes(x=bla, y = estimate, shape = model)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_linerange(aes(x = bla, 
                     ymin = conf.low_95,
                     ymax = conf.high_95), position = position_dodge(width = 1/2), alpha = .7, color = "grey40") +
  geom_linerange(aes(x = bla, 
                     ymin = conf.low_90,
                     ymax = conf.high_90), position = position_dodge(width = 1/2), linewidth = 1, alpha = .7, color = "grey40")   +
  geom_point(aes(x = bla, 
                 y = estimate), position = position_dodge(width = 1/2), color = "black", size = 4) +
  ggtitle("") + ylab ("") + xlab("") +  guides(shape = guide_legend(reverse=TRUE)) + 
  coord_flip()  + #scale_y_continuous(breaks = 0) +
  theme_classic(base_size = 20) + theme( axis.text.y=element_blank(), axis.ticks.y=element_blank()) +
  labs(shape="Measure")  + scale_shape_manual(values = c(56:49)) +  facet_grid(vars(type), vars(fct_rev(term)), scales = "free")



